clear
clear matrix
clear mata
set matsize 4000
set maxvar 10000


*Note: first change working directory replication package folder

run "scripts/programs/program_distlag.do"

/*********************************************************************************/
*Set program parameters
local radius = 20
local angle = 15
local length = 1
/**********************************************************************************/

/************************************** HU Data ***********************************/

use "data/`radius'-by-`length'km_`angle'deg/data_centroid.dta", clear

local controls lnt lnint majorroad tmissing wind_diff
local lags any_res count_res avgval_res val_res habpolytot wsimptot campground wilderness
local costlags tritot roadpcttot

label var roadpcttot "Pct. $< 0.5$ km from road"
label var tmissing "$\Delta T$ missing"
label var lnt "Ln($\Delta T$ + 1)"

local cond `"complexfire != 1 & year > 2011"'

distlag y `controls'  if `cond', ///
	method(almon) params(2,5) lagvars(`lags') costlags(`costlags') ///
	burnt(centroid) tmethod(log) fe(firenum) cluster(firenum) output(table)

gen sample = e(sample)

/********************************* POLICY EXPERIMENTS ***********************************/

qui {

local lead = 1

/************************************* SCENARIO 1 ***********************************/
/********************************* Initial values at 0  ******************************/

mat def results = [0,0,0,.,.]

/* Use precise means
summ avgval_res if sample == 1 & any_res == 1
local init_avgval_res = `r(mean)'
summ count_res if sample == 1 & any_res == 1
local init_count_res = `r(mean)'
local init_val_res = `init_avgval_res'*`init_count_res'
*/

*Approximate for cleaner values
local init_avgval_res = 0.2
local init_count_res = 10
local init_val_res = `init_avgval_res'*`init_count_res'

local P = 2

forvalues p = 0/`P' {
	local exp `exp' _b[wany_res`p']*(`lead'^`p') +
	}
foreach var in avgval_res count_res val_res {
	forvalues p = 0/`P' {
		if `p' == 0 local exp `exp' (
		local exp `exp' _b[w`var'`p']*(`lead'^`p') +
		if `p' == `P' local exp `exp' 0)*`init_`var'' +
		}
	}
noi nlcom `exp' 0
mat v = r(V)
local se = v[1,1]^(1/2)
mat b = r(b)
local b = b[1,1]

mat def results = [results \ `init_count_res',`init_avgval_res',`init_val_res',`b',`se']


forvalues scenario = 2/3 {

	/************************************* SCENARIO 2 ***********************************/
	/******************* Initial values at means within populated cells  ****************/

	if `scenario' == 2 {

		local exp
		/* Use precise means
		qui summ avgval_res if sample == 1 & any_res == 1
		local init_avgval_res = `r(mean)'
		qui summ count_res if sample == 1 & any_res == 1
		local init_count_res = `r(mean)'
		local init_val_res = `init_avgval_res'*`init_count_res'
		*/

		*Approximate for cleaner values
		local init_avgval_res = 0.2
		local init_count_res = 10
		local init_val_res = `init_avgval_res'*`init_count_res'

		mat def results = [results \ `init_count_res',`init_avgval_res',`init_val_res',.,.]

		}

	/************************************* SCENARIO 3 ***********************************/
	/************************ Initial values about 1 SD above  mean  ***********************/

	if `scenario' == 3 {

		local exp
		/*
		local init_avgval_res = 300000
		local init_count_res = 40
		local init_val_res = `init_avgval_res'*`init_count_res'
		*/

		*Approximate for cleaner values
		local init_avgval_res = 0.3
		local init_count_res = 20
		local init_val_res = `init_avgval_res'*`init_count_res'

		mat def results = [results \ `init_count_res',`init_avgval_res',`init_val_res',.,.]

		}

	/************************************* EXPERIMENT 1 ***********************************/
	/************** Change avg. value while holding number of props. constant **************/

	local deltas 2e-1
	local j = 0

	foreach delta in `deltas' {

		local j = `j' + 1
		local comp = `delta'*`init_count_res'

		nlcom (_b[wavgval_res0] + _b[wavgval_res1]*(`lead'^1) + _b[wavgval_res2]*(`lead'^2))*`delta' + (_b[wval_res0] + _b[wval_res1]*(`lead'^1) + _b[wval_res2]*(`lead'^2))*`comp'
		mat v = r(V)
		local se = v[1,1]^(1/2)
		mat b = r(b)
		local b = b[1,1]
		local new_avgval_res = `delta' + `init_avgval_res'
		local new_val_res = `init_val_res' + `comp'

		mat def results = [results \ `init_count_res',`new_avgval_res',`new_val_res',`b',`se']
		}

	/************************************ EXPERIMENT 2 ***********************************/
	/****************** Change count while holding avg. value constant *******************/

	local deltas 10
	local lead = 1

	foreach delta in `deltas' {

		local j = `j' + 1
		local comp = `delta'*`init_avgval_res'

		nlcom (_b[wcount_res0] + _b[wcount_res1]*(`lead'^1) + _b[wcount_res2]*(`lead'^2))*`delta' + (_b[wval_res0] + _b[wval_res1]*(`lead'^1) + _b[wval_res2]*(`lead'^2))*`comp'
		mat v = r(V)
		local se = v[1,1]^(1/2)
		mat b = r(b)
		local b = b[1,1]
		local new_count_res = `delta' + `init_count_res'
		local new_val_res = `init_val_res' + `comp'

		mat def results = [results \ `new_count_res',`init_avgval_res', `new_val_res',`b',`se']

		}

	/************************************** EXPERIMENT 3 ***********************************/
	/************** Change avg. value while holding total value constant **************/

	local deltas 2e-1

	foreach delta in `deltas' {

		local j = `j' + 1
		local comp = ((`init_avgval_res'*`init_count_res')/(`init_avgval_res' + `delta')) - `init_count_res'
		nlcom (_b[wavgval_res0] + _b[wavgval_res1]*(`lead'^1) + _b[wavgval_res2]*(`lead'^2))*`delta' + (_b[wcount_res0] + _b[wcount_res1]*(`lead'^1) + _b[wcount_res2]*(`lead'^2))*`comp'

		mat v = r(V)
		local se = v[1,1]^(1/2)
		mat b = r(b)
		local b = b[1,1]
		local new_avgval_res = `delta' + `init_avgval_res'
		local new_count_res = `init_count_res' + `comp'

		mat def results = [results \ `new_count_res',`new_avgval_res',`init_val_res',`b',`se']

		}

	noi mat list results
	}
}


/************************************* MAKE TABLE ***************************************/

mat colnames results = "count" "avgval" "val" "b" "se"

preserve

svmat results, names(col)
keep count - se
drop if count == .

gen pval = 1-normal(abs(b/se))
gen stars = cond(pval < .01,"**",cond(pval < .05,"*",cond(pval < .1,"+","")))
tostring b, replace format(%9.2gc) force
tostring se, replace format(%9.2gc) force
replace b = "" if b == "."
replace se = "" if se == "."
replace b = b + stars if b != ""
replace se = "(" + se + ")" if se != ""
foreach var in count avgval val {
	tostring `var', replace format(%9.1gc) force
	replace `var' = "{\bf " + `var' + "}" if b == ""
	}
drop pval stars

gen label = ""

replace label = "\T {\bf Scenario I: Initial values}" if _n == 1
replace label = "\T {\bf Scenario II: Initial values}" if _n == 3
replace label = "\T {\bf Scenario III: Initial values}" if _n == 7
replace label = "A. Increase variables to mean within populated cells" if _n == 2
foreach start_val in 1 4 8 {
	if `start_val' != 1 {
		replace label = "A. Increase avg. value while holding no. props. constant" if _n == `start_val' + 0
		replace label = "B. Increase no. props. while holding avg. value constant" if _n == `start_val' + 1
		replace label = "C. Increase avg. value while holding total value constant" if _n == `start_val' + 2
		}
	}
order label

replace se = se + "\\"

local dir "Results/Policy_Expmt/"
capture mkdir `dir'
export delimited using "`dir'table_contents.tex", delim("&") novarnames replace

restore
